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Abstract 



We introduce a geometric method for identifying the coupHng direction between two dynamical systems based on a bivariate ex- 
tension of recurrence network analysis. Global characteristics of the resulting inter- system recurrence networks provide a correct 
discrimination for weakly coupled Rossler oscillators not yet displaying generalised synchronisation. Investigating two real-world 
palaeoclimate time series representing the variability of the Asian monsoon over the last 10,000 years, we observe indications for a 
considerable influence of the Indian summer monsoon on climate in Eastern China rather than vice versa. The proposed approach 
^ can be directly extended to studying K > 2 coupled subsystems. 
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1. Introduction 

Uncovering causal interdependencies from observational data 
is one of the great challenges of nonlinear time series analysis 
ifn . While detecting statistical interrelationships between in- 
teracting systems is comparably simple and can be achieved 
by many linear as well as nonlinear methods such as cross- 
correlation, mutual information fJl, or synchronisation analy- 
sis 1 3 1, the correct identification of coupling direction is still 
faced with numerous practical problems 1 1 1. Existing methods 
for this purpose include linear Granger causality 01 as well as 
nonlinear extensions thereof |5 , 6|, information-theoretic quan- 
tities like transfer entropy |7, 8 | or conditional mutual infor- 
mation 121, and state-space based characteristics such as con- 
ditional probabilities of recurrences (TOlIIIl, to mention only 
a few examples. In general, all methods proposed so far make 
explicit use of dynamical characteristics, and it is still a matter 
of scientific debate which method performs best under specific 
conditions. In turn, the potential geometric consequences of 
directed couplings in the shared phase space of interacting dy- 
namical systems have not been explicitly studied so far. This 
work presents a first attempt to a corresponding structural char- 
acterisation. 

In the last years, several approaches for analysing time se- 
ries from dynamical systems by means of graph-theoretical con- 
cepts have emerged fill [HI El El [El (see lEIEl for a cor- 
responding review). Among these different approaches, net- 
works based on the mutual proximity of state vectors in phase 
space have become increasingly popular, since their structural 
properties are closely related to the geometry of the underlying 
attractor and, hence, the resulting dynamics. 



One particularly important class of proximity networks in 
phase space are recurrence networks (161 [Tvl [Tsl, which are 
especially useful for characterising the structural properties of 
low-dimensional systems. Recall that the property of recur- 
rence (in the sense of Poincare) corresponds to geometrical close- 
ness in phase space without the necessity of trivial temporal 
correlations. This closeness can be characterised in different 
ways, including neighbourhoods of individual states with fixed 
probability mass (/c-nearest neighbour approach) or with a fixed 
phase space volume (determined by a maximum spatial dis- 
tance £ of neighbouring states). In the latter case, the recurrence 
properties obtained from a particular realisation (i.e., a - possi- 
bly multivariate - time series) {xi}f^i representing the relevant 
degrees of freedom of a complex system X are encoded in the 
binary recurrence matrix 



Rij{e) = R (xi.Xjle) = Q{e — d{xi,Xj)), 



(1) 



where <i(-, •) measures some distance (e.g., according to the Eu- 
clidean or maximum norm) in phase space, and 6(-) denotes 
the Heaviside function. The visualisation of this symmetric ma- 
trix is known as the recurrence plot |19|, and the emergence 
of line structures in such recurrence plots has been intensively 
utilised for characterising the dynamical properties of the un- 
derlying time series by estimates of dynamical invariants and 
novel measures of complexity (recurrence quantification anal- 
ysis, RQA) 1201 . 

Recurrence network analysis reinterprets the recurrence struc- 
ture of a time series as the connectivity pattern of an associated 
complex network represented by an undirected simple graph 
|[T6l . Specifically, given the definition in Eq. ([T]) based on e- 
recurrences, we can formally write 
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(where d^j is Kronecker's delta) to obtain the adjacency matrix 
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of the corresponding e-recurrence network (we will abbreviate 
this term by RN in the following). The properties of such net- 
works have been widely studied elsewhere 1 17 . 181 1211 [22ll23l 
[24l|25J, and their practical use as an exploratory tool of time se- 
ries analysis has been demonstrated |[T6l[26l[221. In particular, 
global properties of RNs have been shown to trace dynamical 
transitions in non- stationary time series |tr6l[T8l[26l[23|. A par- 
ticularly relevant, yet straightforward characteristic is the edge 
density 

which equals (for a RN) the recurrence rate RR of the under- 
lying time series. We emphasise that there is a simple mono- 
tonic relationship between p and the recurrence threshold e. 
In the context of RNs (23] |28 1 (but also other spatially em- 
bedded graphs such as climate networks |^3E|), it has been 
argued that a suitable resolution of small-scale network fea- 
tures requires the choice of a low edge density, typically p < 
0.05, confirming earlier suggestions provided for recurrence 
plots 1201 ISTl [321. 

In previous research, the concept of RNs has only been ap- 
plied for the analysis of single (uni- or multivariate) time se- 
ries or for the characterisation of individual dynamical systems, 
whereas the underlying recurrence plot concept has been com- 
plemented by bi- and multivariate extensions |20|. In this work, 
we provide one possible framework for bi- and multivariate RN 
analysis and discuss its potentials for detecting the coupling di- 
rection between two or more complex dynamical systems. Al- 
though our considerations are particularly made for RNs, simi- 
lar conceptual developments should be possible for other types 
of recurrence networks (e.g., /c-nearest neighbour networks). 

This paper is organised as follows: In Section 2, we present 
our methodological approach to a joint treatment of multiple 
time series in a RN framework. The associated characteris- 
tics describing some fundamental structural properties of gen- 
eral coupled networks (sometimes also referred to as intercon- 
nected, interdependent, interacting or network of networks 1 33 1 
in the literature) are reviewed. Subsequently, we outline how 
these properties can be used for characterising possible geo- 
metric signatures of coupling between dynamical systems. The 
proposed framework is then illustrated for two diffusively cou- 
pled Rossler systems as a paradigmatic example in Section 3. 
The effects of different dynamical regimes, the emergence of 
generalised synchronisation, and the lengths of the available 
time series on the detectability of the true coupling direction 
are discussed in some detail. In addition, a real-world example 
from climatology is presented in Section 4. Finally, the main 
conclusions from our work are summarised in Section 5. 



2. Inter-system recurrence networks and coupling detection 

2.7. Cross-recurrence plots 



Let {xi}-^^ and {yj}jJi represent realisations of two dy- 
namical systems X and Y observed at distinct times ti (i = 
1, . . . , Nx) and (j = 1, . . . , Ny), respectively, which do not 
need to be the same for both systems. This is, Xi = X{ti) and 



yj = Y{tj). If both systems share a common phase space (i.e., 
X and Y represent the same physical quantities), we can con- 
sider a distance between individual observations (state vectors) 
of both systems as 



= d'''{X{U),Y{t,)) = \\x,-y,\l (4) 

where || • || is a suitable norm (e.g., Euclidean norm) in phase 
space. Then, the cross-recurrence matrix is 1341 

= e{e^^-d^^{x,^y,)). 



(5) 



Unlike the univariate recurrence matrix R, CR is in general 
not symmetric and not even necessarily a square matrix (i.e., 
Xx 7^ Xy is possible). As in the univariate case, the cross- 
recurrence rate 



RR 



XY 



XY 



(6) 



is an important characteristic of the cross-recurrence matrix CR 
that increases montonically with e^^ . 

The statistical distribution of line structures in the associ- 
ated cross-recurrence plot (CRP) can be utilised for defining 
bivariate measures of complexity to study interrelationships be- 
tween possibly coupled dynamical systems (341 [35l. However, 
this point of view typically requires that the sampling of both 
systems is comparable, since otherwise, line structures are dif- 
ficult to define. 

From a conceptional perspective, we have to emphasise that 
for the formal applicability of CRPs, the two studied systems 
must share the same (reconstructed) phase space, i.e., they need 
to have the same set of dynamically relevant variables [ 20ll32l. 
Although there is no formal proof that this prerequisite pro- 
vides a generally applicable sufficient or even necessary con- 
dition, it provides a reasonable working basis that can be con- 
sidered as being typically met when comparing observations 
with the same physical meaning, such as temperature records 
from two distant observatories or EEG measurements from dif- 
ferent electrodes. In turn, in many real- world problems, one 
is rather interested in studying interrelationships between two 
physically different quantities, for which this requirement is not 
fulfilled. However, in many applications, the resulting concep- 
tual problems have been widely ignored by rendering the re- 
spective phase spaces of two systems quantitatively similar by 
means of suitable monotonic transformations of the data (e.g., 
standardisation to zero mean and unit variance [ ST, '361, quantile 
adjustment, or distinguishing positive and negative interdepen- 
dencies (371 [36l). An alternative solution has been proposed by 
using order patterns, which allows us to directly compare time 
series in terms of their local rank order [ 38 1, but destroys much 
of the quantitative information usually contained in recurrence 
plot-based techniques. Consequently, results based on the in- 
tercomparison of two structurally different variables by means 
of CRPs should always be interpreted and discussed carefully. 

2.2. Inter-system recurrence networks 

In the following, we combine the information on recurrences 
within the individual dynamical systems X and Y (i.e., "intra- 
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system" recurrences, Eq. ([T])) with that on their respective cross- 
recurrences (i.e., proximities between pairs of states of both 
systems, Eq. ([5|) in an inter-system recurrence matrix 



IR(£) 



R (^) 



CR^^(£) R^(£) 



(7) 



with £ being a common threshold for defining the spatial prox- 
imity between the state vectors of each as well as both systems. 
IR can also be derived as the standard recurrence matrix of the 
concatenated observations x = {xi^ . . . ^x^xiVii - - - iVNy)- 
Moreover, since CR^^ = [CR^^]^, IR is a symmetric ma- 
trix. 

We emphasise that the adjacency matrix of two coupled net- 
works (see 1391 ) can be written in complete analogy to Eq. ([7]). 
Specifically, for an inter- system recurrence network (IRN) con- 
structed from two coupled complex systems, we can formally 
define 

A(s) = IR(s)-I^, (8) 

where 1^ is the iV-dimensional identity matrix (A^ = Nx + 
Ny)- In the resulting undirected and unweighted simple (no 
self-loops, multiple edges, or edge and vertex weights) network, 
the vertices represent the state vectors of both systems. Due to 
the block structure of Eq. (|7]), the IRN consists of (i) two (uni- 
partite) RNs of the individual systems, the properties of which 
have been exhaustively studied elsewhere 1161 1171 1181 12T1 [281 , 
and (ii) a bipartite network that includes only edges between 
vertices from different systems. In analogy to CRPs, we call 
the latter bipartite graphs cross-recurrence networks. 

In order to quantitatively study IRNs from a general cou- 
pled networks point of view 1391 (see Section |2.3| ), the two 
subnetworks corresponding to the individual dynamical sys- 
tems need to be clearly distinguished, i.e., there have to be 
only few cross-recurrences in comparison to intra-system recur- 
rences (RR^^ < RR^ ^ RR^). Since this practical require- 
ment is not always met by Eq. (|7]), it is necessary to modify the 
framework introduced above by allowing for different threshold 
distances , and e^^. Equations ^ and ^ \ 



IR(s) 



R^(£^) 

CR^^(£^^) 



1 T 



I thus become 
CR^^(£^^) \ 
R^(.^) ) 



(9) 



with 



and 



A(s) = IR(s)-I^, 



(10) 



respectively. In order to render the two individual RNs quan- 
titatively comparable, it is recommended to choose e-^ and 
independently such that RR^ = RR^ 12311281 . 

Since the IRN framework is based on purely geometric con- 
siderations, the underlying time series do not need to have the 
same sampling. Specifically, IRNs do not make use of any time 
information (which is distinctively different in comparison with 
many other existing methods of time series analysis) and exclu- 
sively require a reasonable attractor coverage in the (original or 
reconstructed) phase space by the sampled data. 



2. 3. Characteristic measures for IRNs 

The topological properties of general coupled networks can 
be quantified by a number of measures which represent a gen- 
eralisation of the canonical measures of complex network the- 
ory |[39l . In the following, we briefly review those character- 
istics which are particularly important for characterising IRNs. 
Subsequently, we discuss how these measures can indicate the 
direction of a coupling between certain dynamical systems. For 
a comprehensive study and detailed background of the underly- 
ing analytical theory, we refer to |23 , 39 1. 

Given a general undirected and unweighted simple graph 
G = (VjE) described by the adjacency matrix A, consider a 
partition of the vertex set V into two disjunct subsets Vx^Vy ^ 
V such that Vx^Vy = V and Vx^Vy = 0. Similarly, the edge 
set E is decomposed into disjunct sets Exx^ Eyy^ Exy ^ E 
with Exx U Eyy U Exy = E and Exx H Eyy = 0, Exx H 
Exy = 0, Eyy fl Exy = such that Gx = {Vx.Exx) 
and Gy = (Vy^Eyy) are the induced subgraphs of the ver- 
tex sets Vx and Vy with respect to the full graph G (Fig. [T]). 
Then Exx contains the (internal) edges within the subgraph or 
subnetwork Gx (likewise Eyy those of Gy), while Exy com- 
prises (cross-) edges interconnecting the subnetworks Gx and 
Gy. Furthermore, we can formally define bipartite subgraphs 
GxY = {Vx U , Exy) containing all vertices of the vertex 
subsets Vx and Vy, and the (cross-) edges between these two 
sets. In the specific case of IRNs, the Gx and Gy correspond 
to the intra-system RNs constructed from the systems X and Y, 
whereas Gxy contains the cross-recurrence structure in terms 
of the sets of cross-edges Exy- In the following, the letters 
v^w^p^q shall denote single vertices (Figs. [T]and|2]). 




Figure 1: A possible representation of two mutually coupled networks. Solid 
lines mark internal edges (of edge sets Exx , Eyy) while dashed lines repre- 
sent cross-edges (elements of Exy) connecting vertices of both subnetworks. 



For a given vertex v G Vx , the cross-degree 

qeVY 



(11) 



gives the number of edges which connect v with any vertex in 
subgraph Gy . Taking the average over all v e Vx and normal- 
ising yields the cross-edge density 



^XY 



NxNy ^ 



(12) 



between the subgraphs Gx and Gy. For an IRN, p^^ equals 
the cross-recurrence rate RR^^ between the systems X and 
Y. 
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The local cross -clustering coefficient 



nXY 



1 



E 



AypApqAqy 



(13) 



estimates the probabiHty that two randomly drawn neighbours 
of vertex v e Vx from subgraph Gy are also neighbours. For 
vertices with a cross-degree of zero or one, we define = 
in order to avoid divergencies. By averaging over all vertices 
V eVx,^Q obtain the global cross-clustering coefficient 



Closely related to is the cross -transitivity 



(14) 



(15) 



Y^veVx;p,qeVY ^vp^vq 

which, as an analog to the canonical network transitivity (401 
|4T]|42l, counts the number of "cross-triangles" over the number 
of "cross-triples". It is important to note that both the cross- 
transitivity and the global cross-clustering coefficient are not 
invariant under the permutation X ^ Y, i.e., 7"^^ 7^ j-yx 
and C^^ ^ C^^ (Fig. [2]), which is in fact the foundation of the 
method presented in the following. 
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Figure 2: Two coupled subnetworks. The graph has global cross-clustering 
coefficients of = 0.5 / C^-^ = and cross-transitivities 7"^^ = 1 / 

r^^ = 0. 



2.4. Network signatures of coupling 

For traditional RNs, the canonical network measure of tran- 
sitivity as well as the local and global versions of the clustering 
coefficient have already been recognised as valuable character- 
istics 1 16][T7l|22l|26]|. Specifically, local clustering coefficient 
and global network transitivity can be interpreted as measures 
for the effective local and global dimension of chaotic attrac- 
tors, respectively |21|. Following and expanding these ideas, 
we now take a look at the "bivariate" versions of these mea- 
sures and how they behave under certain imposed constraints, 
in our case exemplified by a variable coupling between two dy- 
namical systems X and Y. 

By taking into account both the univariate and the bivariate 
recurrence properties of the coupled system, cross-transitivity 
and global cross-clustering coefficient are unique in the way 
in which they act as a filter for information provided not by 
the individual systems' intrinsic processes, but by the way in 



Coupling direction 


Expected relation in network measures 


no coupling 


j-XY r^^q-yX qXY r^^QVX 


X ^Y 


^^YX ^XY ^^YX 


Y ^X 


r^Y ^jYX j^XY ^^YX 


X ^Y 


q-XY ^ q-YX qXY ^ qYX 



Table 1 : Expected qualitative behaviour of IRN measures in different coupling 
situations for systems with comparable properties in the absence of (gener- 
alised) synchronisation. 



which they are coupled. Specifically, if the two systems are 
sufficiently similar with respect to their invariant density, the 
geometric effects of coupling can be understood as follows (Ta- 
ble[T]): 

(i) In the uncoupled case, we can consider T^^, T^^, C^^ 
and C^-^ as randomly arising from the invariant densi- 
ties of X and Y without any additional structural com- 
ponent. Therefore, we can expect ^ q-YX 

C^Y ^qYX 

(ii) Now consider a unidirectional coupling of the type y ex 
f{x — y), where / is a monotonic function of either pos- 
itive or negative sign (for example, in the important case 
of a diffusive coupling, f{x — y) = /j^xy{x — y), with 
fiXY being the coupling strength). The strength and di- 
rectionality of this coupling can be varied. Let Xi and Xj 
be two recurrent states in X. If the coupling direction is 
X ^ Y and the coupling is large enough, we are likely 
to also find a state y^inY, which is (cross-) recurrent 
to both Xi and xj, due to the coupling's diffusive nature 
and thus the tendency to "drag" the trajectory of Y to- 
wards X (Fig. [3]). The resulting "cross-triangle" adds to 
the value of both and C^^ according to their defini- 
tion. On the other hand, "cross-triangles" constituted by 
two recurrent states in Y and one cross-recurrent state in 
X are merely coincidental due to the driver-response-like 
coupling. We would thus expect to see > "J^xy 
C^^ > C^^ in case of a unidirectional coupling X ^Y 
and vice versa for the opposite coupling direction. 

To put it differently: From the viewpoint of the driven 
system Y, the driving system X "contracts" more than 
Y from the perspective of X. Because X is in fact un- 
changed (unidirectional coupling), this effect must be due 
to changes of the spatial distribution of states in Y in- 
duced by the coupling. This can be understood by looking 
at the generic form of the governing dynamical equations, 
y oc f{x — y), i.e., the driven system Y tends to minimise 
its distance to the driver X. 

(iii) Finally, for a symmetric bidirectional coupling, the mu- 
tual effects on both systems are comparable and thus lead 
to IRN measures of the same magnitude. The same ob- 
servation should hold if the subsystems become synchro- 
nised (at least in a generalised sense), e.g., in case of a 
sufficiently strong coupling. 
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The above considerations are still mainly qualitative and 
still need to be supported by a detailed analytical theory of the 
geometric effects of coupling on the trajectories of dynamical 
systems. Specifically, the heuristic explanations presented for 
the (rather wide- spread) special case of diffusive coupling have 
to be further refined, and their extension to other types of cou- 
pling (e.g., impulse-coupling) has to be proven. Both aspects 
provide important directions for further research that are, how- 
ever, beyond the scope of this work. 

3. Example: Coupled chaotic oscillators 

In the following we present a numerical analysis of the ef- 
fect of coupling on IRN measures for a paradigmatic and well- 
studied model system. Specifically, we demonstrate how our 
method performs in the detection of the coupling direction be- 
tween two Rossler oscillators B3l that are diffusively coupled 
via their second component: 

= (1 + v)x^^^ + axx^^^ + /iyx(y^'^ - x^^^) 

= hx^x^^\x^^^ -cx) 
^(i)=-(l-%(2)-^(^) (16) 
^(2) = (1 - ^)^(i) + ay^(2) + /ixy (x^^) - y^^^) 
y^^) =hy^y^^\y^^^ -cy), 

where the parameters ax,y, ^x,y and cx,y define the shape 
(or regime) of the respective attractors of the systems X = 
{x^^\x^'^\x^^^) andF = {y^^\y^'^\y^'^^), provides the sys- 
tems with a frequency mismatch, and jiyx-, I^xy denote the 
respective coupling strengths. Note the diffusive nature of the 
coupling term: The larger the coupling strength, the stronger 
the "dragging" of the driven system's trajectory towards the 
one of the driver. We vary the respective coupling strengths in 
such a way that we get unidirectional (one zero, one nonzero) 
or symmetric bidirectional {jixy = l-^yx) coupling situations. 

For distinct values of the coupling strength, we create en- 
sembles of M = 200 realisations of the considered system 





/ \ 
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• '--.^ - _ - ^ 

■ ■■.:/■■■■ Xj ; 
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Xj 

y oc iixy{x - y) 



Figure 3: Schematic illustration of the influence of a diffusive coupling on the 
IRN structure. In the coupled case (right panel), the trajectory of Y is dragged 
closer to X, while X itself remains unchanged. We thus get an additional 
contribution to C^^ and in form of a "cross-triangle". 



by numerically 1451 integrating Eq. (16) with a step size of 
h = 0.01 for a total time T = 5,000 using random initial 
conditions xo^yo G [0, 1] x [0, 1] x [0, 1]. This ensemble ap- 
proach was chosen to investigate the variability of the network 
measures of interest and ensure the reproducability of the re- 
sults. Of the resulting 500, 000 points on each simulated tra- 
jectory, we discard the first 100, 000 as a possible transient 
phase, and then randomly choose = 1, 500 points as a sam- 
ple of state vectors from the attractor - a technique which is 
widely known as bootstrapping |46|. For both subsystems, we 
then infer M = 200 univariate recurrence matrices (e^) 
and R^(£^) according to Eq. ([T]), each with fixed thresholds 
6:^, £^ such that the recurrence rates RR^ = RR^ = 0.03, 
along with M = 200 cross -recurrence matrices CH^^ (s^^) 
with s^'^ = e^-^ such that the cross-recurrence rate RR^^ = 
0.02. These choices may seem somewhat ad-hoc, but are con- 
sistent with results from previous works (201 ESI EH E2l sug- 
gesting recurrence rates RR^ ^RR^ < 0.05 and considering 
that the coupled networks framework is only meaningful when 
p^^ < ^ (i.e., if the entire network has a pronounced 
modular structure such that it can be truly considered a network 
of networks). 

3.1. Inferring the coupling direction 

We follow the procedure described above for studying two 
different regimes of the Rossler system: the phase-coherent 
regime with 

«x,y = 0.16, bx,Y = 0.1, cx,Y = 8-5, 

and the non-phase coherent "funnel regime" with 

ax,Y = 0.2925, bx,Y = 0.1, cx,y = 8-5. 

In both parameter regimes, a frequency mismatch u = 0.02 
is considered to keep the otherwise identical subsystems from 
completely synchronising at very low coupling strengths, a case 
in which the (sub-)network structures would become fairly iden- 
tical and an identification of the coupling direction by such 
means would be impossible. In case of unidirectional coupling 
(Figs. [4] and [5]), we find that for a large interval of intermediate 
coupling strengths, the coupling direction is identified correctly 
according to our hypothesis: For a coupling direction X ^ Y 
(pxY > 0, Pyx = 0), we notice that the ensemble means and 
standard deviations of and are clearly below those 
of C^^ and T^^ , respectively, and vice versa for F ^ X. In 
general, the results obtained for both measures are qualitatively 
similar. However, although they are conceptually related, both 
quantify distinctively different properties of the networks under 
study, i.e., the observed similarity is no trivial result, but reflects 
the robustness of the geometric modifications to the driven sys- 
tem due to the imposed coupling. 

For small coupling strengths, a discrimination is not possi- 
ble, since the geometric signatures of coupling are too weak to 
be resolved by our network analysis. For stronger couplings, 
the cross-network measures become equal again. We attribute 
this to the occurence of generalised synchronisation, which is 
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0.0 0.2 0.4 0.0 0.2 0.4 0.0 0.2 0.4 



(a) Coupling: X (b) Coupling: Y ^ X (c) Coupling: X ^Y 

Figure 4: (Colour online) Coupling analysis for the Rossler system (Eq. fT6) ) in the funnel regime: Ensemble means and standard deviations (error bars) of the 
cross-network measures C^^, 7"^^ (black) and C^^, T^-^ (red) and the four largest Lyapunov exponents. The Lyapunov spectrum has been estimated using 
the Wolf algorithm | 44l. In the unidirectional cases (a) and (b), the coupling direction is correctly indicated for medium coupling strengths (shaded regions) until 
the onset of generalised synchronisation. In the symmetric bidirectional case (c), the interacting network measures are almost indistinguishable. 
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Figure 6: (Colour online) Differences between the cross-transitivities (top) and 
global cross-clustering coefficients (bottom) for two coupled Rossler systems 
in the funnel regime for X ^ Y (left) and Y ^ X (right) coupling. The well- 
expressed Arnold tongue (white) indicates a parameter region with generalised 
synchronisation. 



indicated by the zero-crossing of the second-largest Lyapunov 
exponent (Figs.|4]and[5]). This observation is consistent with the 
case of bidirectional coupling, where the cross-network mea- 
sures remain equal (Figs.|4jc) andjSj^c)). 

3.2. The onset of synchronisation 

Next, we further investigate the conditions related to the 
emergence of generalised synchronisation and, thus, the break- 
down of our method. For this purpose, we expand our analysis 
into a two-dimensional parameter space, varying not only the 
coupling strength, but also the system's frequency mismatch 
u. The considered discrete {u^ /i) -space consists of 41 x 41 = 
1, 681 points (Fig. [6]) . For each of these points, we generate 
M = 100 IRNs and calculate the mean differences of the net- 
work measures (T^^ — T^^) and (C^^ — C^^), and vice 
versa. 

At low frequency mismatches, the systems synchronise at 
lower coupling strengths, leading to a balancing of the cross- 
network measures (i.e., differences between the measures van- 
ish) and thus a failure of the detection of the coupling direction 
already at relatively low values of the coupling strength. As 
a result. Fig. [6] features white regions, indicating zero differ- 
ence in the measures corresponding to the well-known Arnold 
tongues and revealing the onset of generalised synchronisation 
O. This was anticipated and retrospectively justifies our corre- 
sponding choice of a nonzero detuning as used above. 




^ 1.0 

0.9 



X 



0.8 



■■n-. 0.5 



1 ??l 

6 


, 1 1 

1 


2 00 


5 






N 













60 



1 



10 



5 



N 



RR 
RR 



XY 



X 

0.02 
0.01 



★ 

0.03 
0.02 



0.06 
0.03 



0.10 
0.06 



Figure 7: Ratio of correct identifications of the coupling direction X ^ Y 
by means of cross-transitivity and global cross-clustering for the funnel regime 
of the Rossler system depending on the network size. Higher recurrence rates 
yield better distinctions. In all cases N > 300 is a sufficient lower boundary. 



Moreover, in Fig. [6] positive values indicate parameters in 
which the coupling direction was correctly identified (blue re- 
gions). In turn, negative values correspond to a wrongly de- 
tected coupling direction (red regions in Fig.|6]), which appears 
only for very small coupling strength if the driven system has a 
higher frequency than the driver. This observation suggests that 
in the latter case, the considered detuning has a geometric effect 
on the two attractors that overcompensates the actual signatures 
of (weak) coupling. 

3.3. Limits of the proposed framework 

To explore further properties and possible advantages of our 
method, we now strive to find a lower boundary for the number 
of vertices in the IRN necessary for a correct identification of 
the coupling direction. As an example, we again consider a 
unidirectional X ^ Y coupling between two Rossler oscil- 
lators in the funnel regime with jj^xy =0.1 and a frequency 
mismatch of u = 0.02. In Section |3J1 the chosen number of 
A/" = 1, 500 vertices ensured a good discrimination. In order 
to find a lower bound for N we vary the number of samples 
bootstrapped from the integrated trajectory and repeat our pro- 
cedure again, each time calculating M = 100 realisations of 
the trajectory with random initial conditions for generating our 
IRNs. In order to see whether or not a distinction is possible, we 
consider the number of correct and false identifications (Fig. [7]). 
Moreover, the same procedure was repeated for different values 
of RR^, RR^, and RR^^ to also test the robustness of the 
method with respect to a variation of the recurrence threshold 
(Fig.|7|. 

In all cases examined here, for a perfect discrimination (100 
correct, false cases), no more than roughly N = 300 vertices 
were necessary. This is a remarkable advantage for real- world 
applications, where the amount of available data is naturally 
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much more limited than in the case of model systems. The 
distinction is also robust with respect to variations of the re- 
currence threshold, with higher recurrence rates yielding better 
results for shorter time series. 

4. Real- world example: Palaeoclimate variability 

In order to demonstrate the potential of IRNs for the analy- 
sis of real- world time series, we make use of two proxy records 
of climate variability in the Asian monsoon system during the 
last about 10,000 years (Holocene). The Asian monsoon sys- 
tem is not only one of the most important atmospheric circu- 
lation systems, it also has a strong socio-economic impact be- 
cause it affects a major part of the global population. The insta- 
bility of the Asian monsoon has been correlated with agricul- 
tural and cultural prosperity and political unrest in the past | 47 1. 
Hence, understanding the mechanisms that lead to an inter- 
rupted, weakened, or enhanced Asian monsoon is of paramount 
importance in the on-going discussion on global climate change 
and the monsoon as a tipping element BSl . 

The Asian summer monsoon can be mainly divided into 
the Indian and the East Asian monsoon subsystems (ISM and 
EASM), which transport moisture from distinct sources to the 
continent (Fig. [8]). The study of interrelationships between the 
different monsoon subsystems contributes to a better under- 
standing of the entire monsoon system. Variations in oxygen 
isotope ratios ((5^^0) measured in speleothems provide high- 
resolution proxies for monsoonal activity in the past |[5n[52ll . 
Based on such proxy records, significant spatio-temporal changes 
in the EASM due to different global climate regimes have been 
recently found |53|. However, the ISM is less well understood 
with respect to its long-term evolution, which is partially due 
to the lower number of available high-resolution palaeoclimate 
data from the Indian subcontinent. 

Here, we aim to study the prevalent direction of coupling 
between ISM and EASM since the end of the last glacial pe- 
riod. We use the variability of S^^O measured in stalagmites (a 
proxy for the strength of the summer monsoon) obtained from 
the Dongge cave in Southeast China | 51 1 and the Qunf cave in 
Oman [52 1 (Fig.[8|. The Dongge record D represents Holocene 
monsoon variations in the EASM region, while the Qunf record 
Q contains information on monsoonal climate variability in the 
ISM region (Fig. [9] Tab. [2]). We emphasise that there is no for- 
mal evidence that both data sets actually quantify exactly the 
same climate variable, since the same proxy could have been 
partially influenced by different variables in different ways ac- 
cording to the specific local conditions (e.g., evaporation, drip- 
water rate, CO2 degassing, host rock, etc.). However, in agree- 
ment with the paleoclimate literature |[5T]|52l, we are confident 
that using the same type of proxy data from the same type of 
paleoclimate archive allows considering both records as rep- 
resenting approximately the same integral observable (specif- 
ically, the average annual amount of rainfall, which by itself 
stands as a proxy for the strength of the summer monsoon since 
most of the annual precipitation at both studied sites is in fact 
associated with this season). Moreover, since the global cli- 
mate system is continuous in space, we have both time series 
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Figure 9: Palaeoclimate proxy records of S^^O (see text for a detailed descrip- 
tion) from speleothems sampled from the Dongge (a) and Qunf (b) caves. 

representing the dynamics of different components of the same 
dynamical (physical) system. Consequently, we conjecture that 
the processes recorded in both time series are sufficiently inter- 
related and similar to allow for an application of the proposed 
IRN methodology. 

For the two datasets D and Q we first apply detrending by 
subtracting a 100-yr moving- average from all data (for a de- 
tailed description of the corresponding procedure, see 1271 ). We 
emphasise that this is crucial for the interpretation of all further 
results of our analysis, since we restrict ourselves to variability 
on inter- annual to multi-decadal time scales (i.e., scales below 
100 years). As a second step, we use time-delay embedding of 
the data with an embedding dimension of d = 3 in both cases 
but different time delays of r^) = 3 for the Dongge record and 
tq = 2 for the Qunf record. These parameters have been cho- 
sen according to the results of the false nearest neighbors and 
auto-mutual information methods 1541 . The distinct time delays 
account for the time series' different average sampling rates, 
leading to r ^ 13 years on average for both records to make the 
considered state vectors actually comparable (cf. Tab. [2]). Note 
that for the embedding step, we neglect possible variations of 
the sampling rate with time, as has been done before in suc- 
cessful applications of univariate recurrence network analysis 
to palaeoclimate data 1 16, 27, 26 1. This simplification possibly 
leads to modified positions of individual state vectors in the re- 
constructed phase space. However, since the exact timing of all 
data points is not known and the sampling rate does not change 
too much with time, we argue that within the uncertainty of both 
the recorded observable and the timing of measurements, the 
approximated low-dimensional attractor has a shape that does 
not differ markedly from that one would reconstruct when hav- 
ing "perfect" (equidistant and completely certain) observations. 

In the considered example, the auto-correlations of both 
records decay sufficiently fast (i.e., fall below 1/e within one 
sampling time step on average). This allows us excluding se- 
vere effects due to "false" recurrences corresponding to sub- 
sequent observations, which can occur if the sampling is very 
dense in comparison with the typical time- scale on which the 
auto-correlations decay. In such cases, however, removing the 
corresponding sojourn points from the recurrence matrix by in- 
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Table 2: Time interval T (BP = years before 1950), mean sampling time At, number of observations A^, and corresponding reference of the studied climate proxies. 



Record 


T (yr BP) 


At (yr) 


N 


References 


Dongge D 


-50... 8880 


4.2 


IMA 


Wang et al. (2005) f51J 


Qunf Q 


378. . . 10,300 


7.1 


1,405 


Fleitmann et al. (2003) ^ 



troducing a suitable Theiler window would provide a simple 
and feasible solution 1 17 |. Following these considerations, the 
requirement of a "sufficient attractor coverage" could be turned 
into more explicit conditions (such as time series covering at 
least about 10-20 oscillations of a chaotic system, or a compa- 
rable requirement on the total number of data for general real- 
world time series with a typical preferred time scale). 

Having thus obtained time series in a reconstructed phase 
space that is quantitatively comparable for both records, from 
the embedded state vectors we bootstrap 80% of the data and 
calculate the resulting IRNs and their cross-transitivities and 
global cross-clustering coefficients, with RR^ = RR^ =0.03 
and RR^^ = RR^^ = 0.02 as in Section 3. This procedure 
is carried out M = 10, 000 times to evaluate the robustness of 
our results. The chosen size of the bootstrap samples presents a 
trade-off between the requirements of a sufficiently large num- 
ber of different samples and sufficiently large individual sam- 
ples. We emphasise that because the IRN approach does not 
make any assumptions on the timing of observations, we can 
directly apply this method, including the considered bootstrap- 
ping approach, to the embedded time series disregarding their 
irregular sampling and imperfect chronological control. In fact, 
this is an important advantage of our network-based approach 
in comparison with RQA and similar techniques, which rely on 
line structures in a recurrence plot and, hence, uniform sam- 
pling. 

As a result, we find that the bootstrapped distributions of the 
traditional RN transitivities and global clustering coefficients 
of the two respective records display a considerable overlap 



(Fig. 10). This observation indicates that the two climate sub- 
systems (the dynamics of which have been recorded by the 
considered palaeoclimate archives) are sufficiently similar with 
respect to their long-term dynamics so that both can actually 



be viewed as sharing the same phase space. The latter find- 
ing motivates the further consideration of our IRN approach 
for detecting a possible directional coupling between both sys- 
tems. In fact, the resulting distributions of cross -transitivities 
and global cross-clustering coefficients are clearly separated by 
'j'DQ y jQD qDQ y qQD^ j^^^ general observation per- 
sists when interpolating both records to a common and uniform 
sampling before constructing the IRN (not shown), underlin- 
ing the feasibility and robustness of the considered approach. 
This suggests that the coupling direction was Q ^ D during 
the Holocene. Specifically, under the assumptions that (i) Q 
represents the ISM activity and D the EASM dynamics, (ii) 
auto-correlations do not play a significant role on the consid- 
ered time-scales, and (iii) a possible coupling between both cli- 
mate subsystems is of diffusive nature, we can thus argue that 
the monsoon impact at the location of the Chinese cave was 
probably not only controlled by the EASM but also by the ISM, 
whereas the EASM did not have a comparably strong impact on 
the location of the Oman cave. 

A possible directed influence of one monsoon branch to 
locations usually influenced by the other monsoon branch is 
of considerable importance for understanding the dynamical 
mechanisms of the monsoonal system (e.g., as proposed by 
|55 1). The potential directionality of ISM towards the EASM 
region indicated by our results suggests that the ISM branch is 
probably stronger, at least on average during the considered pe- 
riod of time. This could also imply that the Indian ocean (in 
particular the Indian ocean sea surface temperature, the prin- 
cipal pattern of which is the Indian Ocean Dipole, lOD |36l) 
plays not only an important regional role for the Indian sub- 
continent, but also for the entire Asian monsoon system. The 
EASM is strongly influenced by the El Nino/Southern Oscilla- 
tion (ENSO), whereas the ENSO impact on the monsoon rain- 
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fall in India is known to be more subtle and changing over 
time 157] EH |59l. The detected directionality adds a new as- 
pect to this relationship and might help in further improving 
our understanding of the Asian monsoon system as a whole. 
We particularly emphasise that IRNs can provide important in- 
formation for constructing and interpreting networks of palaeo- 
climate records [60| to infer dominating spatial patterns of 
climate dynamics in the Earth's history. 




1.4 0.45 0.5 0.55 0.6 0.65 0.7 



c 


D 




c° 






1 

1 


1 


1 



0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 
C 



Figure 10: (Colour online) Normalised distributions of the IRN measures cross- 
transitivity (top) and global cross-clustering coefficient (bottom) inferred from 
the Dongge and Qunf data by bootstrapping 80% of the data in M = 10, 000 
independent realisations. The distributions of the inter-system measures are 
clearly separated with {Q^Q ^q^DQ ^ qQD q-QD^^^ indicating a possible 
coupling direction Q ^ D. 



5. Conclusions 

We have proposed a new complex network-based approach 
for inferring coupling directions from bivariate time series. The 
recently introduced concepts of recurrence networks ifTSl and 
coupled network analysis 1 39] have been combined to the new 
analytical tool of inter-system recurrence networks (IRNs). Cross- 
transitivity and global cross-clustering coefficient of IRNs can 
be used to quantitatively derive the coupling direction based 
on purely geometric considerations. Although developing a 
corresponding generally applicable theoretical understanding, 
as well as a systematic comparison with other existing tech- 
niques for detecting the coupling direction between two mea- 
sured signals, will be subject of future research, our results 
suggest that the proposed framework can be applied to many 
situations where the exact functional form of coupling is not 
necessarily known. 

We have demonstrated the potential of the proposed frame- 
work using a prototypical model system and discussed its ro- 
bustness and sensitivity with respect to data length and recur- 
rence threshold. We find that the IRN method performs well 
even for short time series (TV ^ 300). Another important ad- 
vantage of our approach is that it relies on the spatial structure 
of the coupled attractors in phase space and does not require any 
time-information. Hence, problems of strong auto-correlations 
often present in real-world systems are avoided by construc- 
tion. In turn, our method is not able to detect possible coupling 
delays unless being applied separately for different, mutually 
shifted time slices. 

The proposed approach can be directly generalised to study- 
ing mutual couplings among a set of i^T > 2 subsystems, which 
leads to IRNs the adjacency matrices of which have a block 
structure composed of K traditional (unipartite) recurrence net- 
works and K{K — l)/2 bipartite cross-recurrence networks. 
Note, however, that the study of couplings for i^T > 2 is much 
more challenging than for i^T = 2, since in this case both direct 
and indirect interactions may be present and need to be distin- 
guished. It will be subject of future research to investigate if 
and how the proposed geometric framework based on multi- 
variate extensions of recurrence networks is able to cope with 
this problem. 

As a real-world example, we have studied the possible di- 
rectionality of coupling between the Indian Summer Monsoon 
(ISM) and the East Asian Summer Monsoon (EASM) during 
the last about 10,000 years. Our method allows its direct appli- 
cation to the data without interpolation to a common time axis. 
From our analysis we find indications that the ISM is probably 
influencing the monsoon strength in Eastern China (located in 
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the EASM influence region) on multi-decadal time scales be- 
tween about 10 years (data resolution) and 100 years (scale of 
detrending), which supports previous findings based on mod- 
elling |55|. However, further research is necessary to unam- 
biguously verify this result by means of complementary meth- 
ods of bivariate time series analysis. 

As a final remark, it is worth mentioning that there are other, 
conceptually different multivariate generalisations of recurrence 
plots, such as joint recurrence plots 116111 and multivariate recur- 
rence and cross-recurrence plots 1621 . In principle, all these 
concepts could be re-interpreted from a complex network point 
of view as well. However, due to their different construction 
mechanisms, they provide different types of information on the 
investigated systems in comparison with the approach presented 
in this work. A detailed corresponding study will be subject of 
future work. 
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